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(57) ABSTRACT 

.An image-based phase retrieval technique has been devel- 
oped that can be used on board a space based iterative trans- 
formation system. Image-based wavefront sensing is compu- 
tationally demanding due to the floating-point nature of the 
process. The discrete Fourier transform (DFT) calculation is 
presented in “diagonal” fonn. By diagonal we mean that a 
transformation of basis is introduced by an application of the 
similarity transform of linear algebra. The current method 
exploits the diagonal structure of the DFT in a special way, 
particularly when parts of the calculation do not have to be 
repeated at each iteration to converge to an acceptable solu- 
tion in order to focus an image. 
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DISCRETE FOURIER TRANSFORM IN A 
COMPLEX VECTOR SPACE 

PRIORITY CLAIM 

5 

This application claims priority under 35 U.S.C. §119 to 
U.S. Provisional Application Ser. No. 61/311,909 entitled 
“DISCRETE FOURIER TRANSFORM IN A COMPLEX 
VECTOR SPACE” filed on Mar. 9, 2010, the entire contents 
of which are hereby incorporated by reference. 1° 

ORIGIN OF THE INVENTION 

The invention described herein was made by an employee 
of the United States Government, and may be manufactured 15 
and used by or for the Government for governmental pur- 
poses without the payment of any royalties thereon or there- 
fore. 

BACKGROUND OF THE INVENTION 20 

1. Field of the Invention 

Hiis application relates to Discrete Fourier Transforms, 
and in particular, to Discrete Fourier Transforms in phase 
retrieval for use in waveform analysis. 25 

2. Background 

Phase retrieval is an image-based wavefront sensing 
method that utilizes point-source images (or other known 
objects) to recover optical phase information. The most 
famous application of phase retrieval was the diagnosis of the 30 
Hubble Space Telescope mirror edge defect discovered soon 
after the launch of Hubble and subsequent correction using 
Corrective Optics Space Telescope Axial Replacement 
(COSTAR.) It may be ironic that a phase retrieval wavefront- 
sensing method lies at the very heart of the commissioning 35 
process for Hubble’s successor. The earliest suggestion of 
using phase retrieval as a wavefront sensing method for 
JWST was in 1989, nearly a year before the Hubble launch 
and deployment. Details documenting the Hubble Space 
Telescope phase retrieval analysis are discussed in the litera- 40 
ture. 

A number of image-based phase retrieval techniques have 
been developed that can be classified into two general catego- 
ries, the (a) iterative-transform and (b) parametric methods. 
Modifications to the original iterative-transform approach 45 
have been based on the introduction of a defocus diversity 
function or on the input-output method. Various implemen- 
tations of the focus-diverse iterative-transform method have 
appeared which deviate slightly by utilizing a single wave- 
length or by varying the placement and number of defocused 50 
image planes. Modifications to the parametric approach 
include minimizing alternative merit functions as well as 
implementing a variety of nonlinear optimization methods 
such as Levenburg-Marquardt, simplex, and quasi-Newton 
techniques. The concept behind an optical diversity function 55 
is to modulate a point source image in a controlled way; in 
principle, any known aberration can serve as a diversity func- 
tion, but defocus is often the simplest to implement and exhib- 
its no angular dependence as a function of the pupil coordi- 
nates 60 

The discrete Fourier transform (DFT) calculation is pre- 
sented in “diagonal” form. By diagonal we mean that a trans- 
formation of basis is introduced by an application of the 
similarity transform of linear algebra. 

What is needed is to find a more efficient implementation 65 
of the DFT for applications using iterative transform meth- 
ods, particularly phase retrieval. 


2 

BRIEF SUMMARY 

The image-based wavefront sensing method of and appa- 
ratus for focusing by determining phase difference between 
received images, removing that phase difference by calculat- 
ing a DFT for each signal and iteratively multiplying the 
resultant until a convergence is achieved comprises the steps 
of accepting point-source stellar objects to recover embedded 
optical phase information, receiving data from a sensor array 
wherein the data includes embedded phase information, cal- 
culating a similarity transform and pre-transforming the data 
once at the beginning of the calculation. 

The image-based wavefront sensing method of and appa- 
ratus for further includes iteratively processing the data using 
arbitrary frequency spacing in a diagonal fashion until a 
wavefront convergence is achieved to determine a phase dif- 
ference by implementing the DFT in a 1 dimensional linear 
array using just N operations, where 0<N<Nxlog(N)<oo, 
wherein Nxlog(N) is the number of operations of an FFT 
counterpart, back-transforming the data at the end of the 
calculation in a single step and focusing the received images 
based on the determined phase. 

The general problem of diagonalizing the DFT matrix 
(DFT eigen-problem) has been around for years. Apparently 
some aspects of the problem are related to topics in number 
theory that were even solved by Gauss. Such considerations 
have not led to more efficient implementations of the DFT 
over the FFT implementation, but we note that for certain 
applications, e.g., phase retrieval, that it is possible to exploit 
the diagonal structure of the DFT in a special way, particu- 
larly when parts of the calculation do not have to be repeated 
at each iteration to converge to an acceptable solution. 

BRIEF DESCRIPTION OF THE SEVERAL 
VIEWS OF THE DRAWINGS 

FIG. 1 is a functional block diagram of a system which 
contains a discrete Fourier transform processor in a complex 
vector space, according to an example embodiment; 

FIG. 2 is a flowchart representing the methodology of a 
discrete Fourier transform processor in a complex vector 
space, according to an example embodiment; and 

FIG. 3 is another flowchart representing the methodology 
of a discrete Fourier transform processor in a complex vector 
space, according to an example embodiment. 

DETAILED DESCRIPTION 

Detailed illustrative embodiments are disclosed herein. 
However, specific structural and functional details disclosed 
herein are merely representative for purposes of describing 
example embodiments. Example embodiments may, how- 
ever, be embodied in many alternate forms and should not be 
construed as limited to only the embodiments set forth herein. 

Accordingly, while example embodiments are capable of 
various modifications and alternative forms, embodiments 
thereof are shown by way of example in the drawings and will 
herein be described in detail. It should be understood, how- 
ever, that there is no intent to limit example embodiments to 
the particular forms disclosed, but to the contrary, example 
embodiments are to cover all modifications, equivalents, and 
alternatives falling within the scope of example embodi- 
ments. Like numbers refer to like elements throughout the 
description of the figures. 

It will be understood that, although the terms first, second, 
etc. may be used herein to describe various elements, these 
elements should not be limited by these terms. These terms 
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are only used to distinguish one element from another. For 
example, a first element could be termed a second element, 
and, similarly, a second element could be termed a first ele- 
ment, without departing from the scope of example embodi- 
ments. As used herein, the term “and/or” includes any and all 5 
combinations of one or more of the associated listed items. 

The terminology used herein is for the purpose of describ- 
ing particular embodiments only and is not intended to be 
limiting of example embodiments. As used herein, the singu- 
lar forms “a”, “an” and “the” are intended to include the plural to 
forms as well, unless the context clearly indicates otherwise. 

It will be further understood that the terms “comprises”, 
“comprising,”, “includes” and/or “including”, when used 
herein, specify the presence of stated features, integers, steps, 
operations, elements, and/or components, but do not preclude 15 
the presence or addition of one or more other features, inte- 
gers, steps, operations, elements, components, and/or groups 
thereof. 

Hereinafter, example embodiments will be described with 
reference to the attached drawings. 20 

Example embodiments of the present invention include a 
Discrete Fourier Transform (DFT) in phase retrieval for use in 
waveform analysis in a space based application. 

One such space based application under consideration to 
use DFT phase retrieval methods includes the James Webb 25 
Space Telescope (JWST), one of NASA’s great observatories 
of the Origins program and is scheduled to launch later this 
decade. The JWST design includes using advanced optical 
technology and 0-30 grade beryllium mirrors. The primary 
mirror is a 6.5 meter segmented primary using 1 8 segments in 30 
a hexagonal array. The optical design and manufacturing 
tolerances are specified for diffraction-limited performance 
at 2 pm, and the telescope is designed for operation over a 
temperature range of 30-60 K with an areal density of less 
than 15 kg/m 2 . As a result, JWST commissioning and peri- 35 
odic optical maintenance incorporate san image-based wave- 
front sensing and control (WFSC) system to align the mirror 
segments, minimize the effects of figure error, and position 
the secondary mirror of the three-mirror anastigmat design. 
The wavefront sensing method specified for JWST is image- 40 
based in the sense that point-source stellar images are col- 
lected and processed to recover optical phase information. 
The primary camera for this function is the JWST Near Infra- 
red Camera (NIRCam), although it has been shown that the 
Near Infrared Spectrograph (NIRSpec) is suitable for this 45 
function as well. 

FIG. 1 shows a block diagram of a satellite or platform in 
which an embodiment of the instant invention is contained. 
The sensor array 130 is connected both physically and elec- 
tronically to the platform housing 120 which is in turn like- 50 
wise connected the processor 110 which contains the primary 
embodiments of the instant invention. The platform housing 
includes all the necessary parts (not shown) of a satellite such 
as fuel, electrical power, communications etc. Image-based 
wavefront sensing is computationally demanding due to the 55 
floating-point nature of the process. With certain data con- 
figurations and large array sizes, and particularly for under- 
sampled data sets, the calculations can take tens of minutes or 
even hours to reach full conveigence for the optical wave- 
front. In addition, processors are generally susceptible to 60 
opto-mechanical instabilities induced by environmental fac- 
tors such as jitter and turbulence. These factors not only limit 
fine-phasing performance but also the stability of the phasing 
results. In such cases, an obvious mitigation strategy is to 
increase the computational bandwidth of the WFSC system, 65 
i.e., decrease the time required for an entire closed-loop 
WFSC cycle. 


To improve conveigence times without sacrificing accu- 
racy and to mitigate these environmental factors, we have 
developed special-purpose floating-point intensive comput- 
ing hardware. These processors communicate using a series 
of multi-processor computing architectures specifically tai- 
lored to the JWST fine-phasing methodology. The computing 
system is based on a cluster of 64 DSPs (digital signal pro- 
cessors) operating in parallel under a shared memory archi- 
tecture. The system is portable and callable from a laptop 
computer using a typical Ethernet connection, reaching a 
sustained floating-point performance of approximately 100 
Giga-Flops with 64 DSPs for development, testing and hard- 
ware redundancy. 

An embodiment of the instant invention is illustrated by 
some practical aspects of calculating in a diagonal DFT basis. 
Aside from the application to iterative Fourier methods, cal- 
culating in the diagonal DFT basis also has an interesting 
geometrical interpretation as solving for the Fourier trans- 
form pair, (F„, f„), that are co-linear or “parallel” with one 
another. Additional discussion on this interpretation is given 
below. 

As an overview we begin by writing down the matrix form 
of the DFT. Next the hybrid-index notation of Schouten is 
introduced as an efficient notation for performing the DFT 
calculations in a complex vector space. The similarity trans- 
formation is presented in detail for a low-dimensional DFT 
example (4x4 array size) and then applied to illustrate how the 
DFT can be carried out in a diagonal basis using the hybrid- 
index notation. The main results show that by “pre-transform- 
ing” the data once at the beginning of the calculation, and then 
“back-transforming” the data at the end of the calculation, the 
DFT can potentially be implemented as a 1 -d vector product 
using just N operations, compared to the Nxlog(N) operations 
needed by the FFT counterpart. The “pre-” and “back-” trans- 
forms are implemented using the DFT similarity transform 
and we further note that the calculation can be implemented 
using arbitrary frequency spacing. 

DFT as a Unitary Transform 

The DFT calculation can by implemented in matrix form: 


N-l 

F(v m ) S DFT(v„) = fO„)V m „, 

n = 0 


( 1 ) 


where V mn is a non-singular MxN Vandermonde array: 

y =e -i 27 tv rt ^nAx)— e -i 2 Tximn)IM =: ^ e -i 2 n/Mynn = ^ > ^mn^ 

with components given by 


Vmn 


1 1 
1 OJm 

i <4 # 

M M 
1 oj m 
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CO 


2 
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CO 


2-2 
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L 1 

/ / i N ~ l 

L OJ M 

i , 

L Cx) M 

O M 

, CM- lXA'-l) 

L OJ M 


(3) 


Next we define W mw as 

w„„= v mn MM={ e - a * IM Y' -Mm, 


(4) 


and then it is straightforward to show that W m „ is a unitary 
matrix: 


w-W<=l, (5) 

using boldface matrix notation where I is the identity matrix 
and the “f ” symbol denotes the complex-conjugate transpose 
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or “adjoint.” We co mm ent that the adjoint operation will be 
presented more generally in the next Section. Next letting 

( 6 ) 

Equation (1) takes the form of a unitary transformation: 5 

N-l (7) 

F(v m ) = F„ = £ W™/„, 

n=0 

with the inverse transformation given by (generally but not to 
always we take M=N): 

M-l (8) 

fn = 2 WnmfFn. 
m= 0 

15 

Given (5) and (7), the DFT can be analyzed using the 
transformation theory of complex vector spaces. The hybrid- 
index notation of Schouten supplies a natural framework for 
this approach and is presented in the next Section. 

Hybrid Index Notation 20 

As illustrated above, the DFT is a unitary transformation in 
a complex vector space. Therefore, the calculations can be 
performed elegantly and efficiently using the tensor-transfor- 
mation formalism of Schouten. In particular, Schouten’s 
hybrid-index notation provides a natural framework for the 25 
DFT calculations. The goal in using this approach is to give a 
more general representation of the calculations by identifying 
components with quantities of a geometric manifold. Often, 
but not always, such identification can lead to alternative or 
more efficient calculation strategies as discussed later. There- 30 
fore, the goal in the next Section is to sort out the details of the 
hybrid-index notation for the DFT, and thus identify the DFT 
with quantities in a geometric manifold. In addition, by 
adopting a notation that makes the transformation rules of a 
quantity explicit, the meaning of the calculations can often be 35 
clarified and the notation serves as an efficient bookkeeping 
device for checking consistency throughout the calculations. 

DFT in Hybrid-Index Notation 

The DFT is defined here as a unitary transform between 
conjugate Fourier domains. The 1-d DFT transform is then 40 
expressed using Schouten’s hybrid-index notation as 


A point transformation, on the other hand, is denoted by 

ir=p.„ M i ( 12 ) 

and is the form adopted for the DFT in (9). Another important 
distinction between a transformations of basis and a point 
transformation is that the indices for the point transformation 
are ordered and the positions are marked using over- or under- 
dots as in P.„ m (12). The main take-away is that for coordinate 
transformations, the kernel symbol, u, remains the same as in 
(11). But for point transformations, the kernel symbol 
changes (as in (12)) but the indices remain the same. So 
initially at least in this context, the DFT transformation in (9) 
is defined as a point transformation, but later we find that 
when calculating the DFT in a diagonal basis, it is natural to 
reconsider (9) using a common kernel symbol for the Fourier 
transform pair, specifically, (F m , f” )-»(f-, f ). We note also 
that it is possible to mix the two transformations (basis and 
point) as is required when considering the allowable transfor- 
mations of basis in digital signal processing. These details 
will be discussed in more detail later. 

Continuing with an explanation of the formalism, the com- 
plex conjugate of (9) is given by 

F"=n";f=f- (13) 

using the over-bar to denote the complex conjugation of com- 
ponents and kernel symbols in accordance with Schouten’s 
hybrid-index notation and noting for consistency that 

F 5 ”=F"=>F"=E' m =F". (14) 

The conjugation of indices and kernel objects is fundamen- 
tal to the study of transformations in a complex vector space, 
and as illustrated below, is perfectly suited for calculations of 
the DFT transform. 

We choose an initial convention that the data values are 
labeled by contravariant components, f", and then by com- 
plex conjugation of the kernel symbol and indices in various 
combinations, several other types of components are repre- 
sented in Schouten’s notation by the “raising” and “lowering” 
of indices using the “fundamental” or metric tensor, a a : 


F m =W m .J", 


(9) 


fata. 


(15) 


using the Einstein sum convention defined over repeated indi ■ 
ces in an upper/lower combination: 


where the latter two follow by complex conjugation of the 
45 first two. Specifically, we mean that 


2 W m „f„ = 


(10) 


The lower/upper indices denote the covariant and contra- 
variant transformation rules — the upper indices are not expo- 
nents in this notation. Exponents of a quantity will be denoted 
using parentheses, or brackets, respectively, for either a trans- 
formation of coordinates, or a “point transformation.” Allow- 
able transformations of coordinates for this application will 
be discussed later. The distinction between the transforma- 
tions of coordinates and point transformations are fundamen- 
tal to Schouten’s kernel-index method. To elaborate briefly, 
coordinate transformations are denoted by a transformation 
of basis, say, m—»m\ where the new basis set is denoted by m'. 
The components in the new basis are thus 

u m '=A m m 'u m , (11) 

where A m M is the transformation matrix and the U™ are the 
vector components in the new basis set (u is arbitrary). 


r=>Ar‘>sJ"3"=>J m =‘>sJ", (16) 

As mentioned above, co- and contravariant transformation 
_ n rules can be switched back and forth using the metric or 
fundamental tensor, a B , to raise or lower indices as appropri- 
ate. The metric tensor also defines how the differential ele- 
ment of distance is calculated in the manifold as well as the 
inner product between two vectors, say, u A and v 7 : 

(u\v) u ; \f tlvi. (17) 

pointing out that in a linear complex vector space, the “bra- 
ket” notation of Dirac is easily incorporated in this approach 
as illustrated on the LHS of (17). Illustrating the use of the 
60 metric for raising and lowering indices we have: 

and/^/i, (18) 

and then substituting these expressions back into (9) we see 
that an alternative form of the DFT is obtained: 

65 

F™ = W'r=(aijrHa"%=a- ln a^W™%= 


( 19 ) 
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Adjoint in a Geometric Manifold 

In general, the adjoint operation applied to the W.,,”, 
denoted by is defined through a two-step 

operation using both the metric tensor and complex conjuga- 
tion: 5 

(a) each index is raised or lowered as appropriate using the 
metric: 

( 20 ) 

(b) the result in (a) is then complex conjugated to give io 

(W.„ m ) f = a * k a ' m ^*1 = W'- m =w„ m , (21) 

and finally we write 

W. n m =={W. n m )1='W„ m , (22) 

For comparison, in matrix boldface notation: 

(W) x =W r =aWcr l = Vr l , (23) 

where “T” defines the usual transpose operation. We note 
especially that the example in (23) illustrates a degree of 20 
non-generality in the way tire adjoint is usually calculated in 
matrix notation. Specifically, defining the adjoint using com- 
ponents in a geometric manifold (using steps (a) and (b) 
above) can give different results than simply complex conju- 
gating and then performing the transpose operation using the 25 
“T” as in (23). In special coordinates, a- k =b^ k , these two are 
identical, but in general, that is not always the case for arbi- 
trary a at- 

Since the W.„ m are unitary, the inverse, 


r 
Tr ■«’ 

is defined using the adjoint and is expressed (from (22)): 35 


_i (24) 

w m „ = w m n = 

40 

such that 

W. k "W=b.l (25) 

We further comment that Schouten did not initially use an 
over-bar on the kernel symbol as in we have in (25). Consid- 45 
ering the W.„ m using two upper indices, W"'“ (see (18) and 
(19)) the corresponding adjoint is similarly calculated: 


8 

scaling and phase transformations. The details for these vari- 
ous transformations in this notation will be discussed in a later 
set of notes. 

Matrix Multiplication Using Index Notation 
When calculating with components, index order is impor- 
tant for proper matrix multiplication and thus we adopt the 
following convention for a rank-2 quantity: 


A™, (29) 

1 


and the inverse of A.„ m is notated as, 


-1 



such that 


A m „A" k =6 m k 

where b. k m is the Kronecker delta. To illustrate, the matrix 
multiplication of two matrices A and B is notated as: 

c.r=A,rB. k \ ( 30 ) 

and to be specific we calculate (30) using Mathematica with 
2x2 arrays: 

7Me[S\\m[(A[[um,n]]B[[n,lk]]),{n,\,NN}],{um,\,2}, 

{/A; 1,2}]] (31) 

where um denotes “upper-m” and lk denotes “lower-k;” the 
upper/lower dummy variable, n, is left un-notated in Math- 
ematica with regard to its position. The result is of course 
(30): 


Cl = AZE? k =AB = 


a\ b\ + a\b\ a\ b\ 4 - a\b\ 
a\b\ 4 - a\b\ a\ b\ + a\b\ 


(32) 


where AB is the conventional matrix product using matrix 
boldface notation. For comparison, transposing both the m-n 
and n-k index positions in (30) reverses the order and gives the 
BA matrix product: 


(W m7I ) t = a mt a lrl W tl -W mn = = Wj. 


(26) 


The inverse DFT transformation to (9) is thus defined: 

W. m "F" = ww. k m /6 k "f=f, (27) 

as expected, or when using the form of (19) the inverse DFT 
is expressed 

( 28 ) 

A few additional comments from this Section: in tensor 
transformation theory, there is geometric meaning attributed 
to these differing but complementary transformation rules. 
Specific transformations are defined by an allowable Group 
of transformations, G a , as discussed further in Schouten. The 
corresponding allowable transformations of basis for the 
DFT application (digital signal processing) can include frac- 
tional and integer shifts of the data origin in addition to 


50 


Cl = A™B? =BA = 


b\a\ + b l 2 a k 
b\a\ + b\a\ 


b\ a\ 4 - b\ a\ 
b\ a l 2 4 - b\ a 2 . 


(33) 


55 and the corresponding Mathematica code for (33) is 

Table[Sum[(A/y«, um]]B[[lk,n\ ] ), {n, 1 JSIN}], {um, 1 ,2 }, 

{lk, 1,2}] (34) 

Metric Tensor Deriving from Unitary Invariance 
60 In summary of Schouten’ s notation, various geometric 
quantities are defined that are distinguished by their transfor- 
mation rules using index position. A group of allowable trans- 
formations, G a , is introduced and then several invariant quan- 
tities may be constructed based on invariance of a given 
65 symmetry group. One such quantity is the metric tensor, a B , 
defining “distance” in the manifold and also defining the 
inner product between any two vectors . Therefore, one goal in 
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these notes is to solve for the DFT metric tensor, a B , by 
postulating nothing more than unitary invariance under the 
action of the point transformation, W.„ m , as illustrated below. 
Specifically, given that the W.„ m are defined as a unitary 
transformation, we introduce a metric tensor defined by 5 
invariance under the action of the W.„ m transformation such 
that 

w.rwfa- m =a- kh (35) 

or in matrix form: io 

WaW~ l =a. (36) 

Quantities that are invariant under the action of the unitary 
group are hermitian and thus the a B satisfy 

a B =(aB)'=a ]t , (37) 15 

The introduction of a hermitian metric a B into the E„ geomet- 
ric manifold defines the manifold as a U„ which differs from 
the earlier notation used by Schouten where this manifold was 
labeled as the R M . 20 

General form of Parseval’s Theorem 

We can show that Parseval’s theorem is really implied by 
unitary invariance, and not necessarily “energy conservation” 
as it is commonly referred to in the engineering literature. To 
begin, we consider the magnitude squared of F"'-» I F 1 2 , giving , 5 

(F\F)=\F\ 2 =a li J'”F'. (38) 

Substituting the DFT transformation rule from (9) into (38) 
we have 


(39) 30 

but since the a B are invariant under this unitary transformation 
(see (35)) then the RFtS of (39) simplifies and we get 

IFI 2 =ai//. (40) 

Therefore, by comparison with (38), Parseval’s theorem is 35 
seen to hold for any solution of the a s as implied by unitary 
invariance (35), and thus in general: 

a T J i F'=a l ff. (41) 

The common form of ParsevaTs theorem is expressed in 40 
vector notation as 


\F\ 2 =\f 2 . (42) 

To give some perspective on these results we look forward 
a bit and comment that the standard DFT approach adopts a 45 
special solution (coordinate system) for the a^, given by the 
Kronecker delta (as shown later): 

(43) 

and so in index notation (42) takes the familiar “Euclidean” 50 
form: 

(44) 

or as it is commonly written: 

\F l \ 2 +\F 2 \ 2 +L+\F N \ 2 =\f l \ 2 +\f 2 \ 2 +L+\f N \ 2 . (45) 55 

The point is that for a general a B , ParsevaTs theorem (41) 
will not necessarily take the Euclidean form as demonstrated 
later below (see (67)), although unitary invariance implies 
that (41) still holds. 60 

A couple of additional co mm ents: as discussed later the 6 
s is a solution for the a B and corresponds to a positive definite 
signature solution. The approach taken here is motivated by 
that fact that in general, the eigen-problem requires the use of 
a general metric, and by no means does this necessarily cor- 65 
respond to the Kronecker delta in a unitary complex vector 
space, U„. But even so, we later adopt the Kronecker delta 


10 

solution to make progress on the main point of calculating the 
DFT in a diagonal basis, and for easy comparison with pre- 
viously published results. 

Metric Solutions 

By postulating unitary invariance of the metric we can 
solve for the a B using (35) and (37). The W. k m are known 
through identification with (4): 


w* = -7=(e" 
VM 



(46) 


As discussed above, the inverse transformation of (46) is 
defined through the adjoint operation in (20). But at this point 
in the discussion we have not yet solved for the metric aj„ 
which is needed to form the adjoint. Therefore, we first solve 
for the 


-i 


W'i 


such that 


KW" k =S%, 


(47) 


and the 


are given by: 


1 ( c i2n (MVMj _ J_( 1 1 

j~m VTU -l, 


(48) 


Continuing to the general 2x2 metric, a B is defined as 


a u - 


Z 4i 

Z T2 

fxn+iyn 

*12 + iy n ' 

Mi 

Z 12 , 

v *21 + iy 21 

*22 + iy22 


(49) 


where the complex notation is dropped on indices of real 
components. The condition for unitary invariance (35) then 
gives the following simultaneous system of 4 equations and 8 
unknowns: 


Z Ji 

+ J T2 +%i +JJ2 

Z T1 ~ Z T2 + % “ %2 ' 

Z Tl 

Z T2 ' 

Z Tl 

K? 1 

1 

tV 

1 

+ 

Z I1 “ Z T2 - % + %2 , 

Z 21 

Z 12 > 


The system (50) is under-determined and can therefore 
have many possible solutions but will be further constrained 
by the hermitian condition (37). Solving (50) eliminates 4 of 
the unknowns (z Xl and z T2 ) in terms of the other 4 unknowns 
(zj| and z 22 ) showing that a B has the general form: 
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+ Zj2 Zji (51) 

Z 21 Z 22 , 


Next applying the hermitian (symmetry) condition (35) to 
(51) gives the condition that the 2x2 a# must be real valued: 


2^21 + yi2 
J21 


^21 

^22 



=> Im(a f/ ) = 0, 


( 52 ) io 



1 * 1 

(i 

0) 


1 * ( 

' 1 1 


^ kl 

|*21-° = | 

io 

J 

’ ^ kl 

hi- 1 , - 

2 1 "I 

) 


*22-* 


*22- _1 








0 

1 ' 


1 * 1 

f 2 

1) 



"2 

kl 

II 

—i o 

h n 

it 

oJ 

l : ?n 

1*21 -“2 
*22-! 

1 

<~2 

1 


( 55 ) 


but as expected there are multiple solutions for the a B in terms 
of the 2 final parameters x 21 and x 22 . The general solution for 
the 2x2 case is thus reduced from 8 to 2 unspecified real 
components and is given by 


2^21 4 •'■22 *21 1 ( 53 ) 

*21 *22 > 

Asacheck, substituting (53) back into (35)verifies that tliis 
s.jj is invariant under the action of the W k m for arbitrary (but 
real) numerical values of x 21 and x 22 in (53). 

We have shown that an infinite number of solutions exist 
for the 2x2 metric, a B , as implied by unitary invariance in 
terms of two real parameters. We note that one of the simplest 
and most important solutions is a special case of (53) when 
x 21 =0 and x 22 =l giving the Kronecker delta as mentioned 
earlier: 


Given some specific solutions for the a B in (55), we examine 
various other quantities that are derived in Schouten’s nota- 
tion. For instance, consider the second solution in (55) that is 
a multiple of the W.„ m : 



and then for generality define a Complex data vector: 

/=(z 1 .z 2 M* 1 +<> 1 .* 2 +‘> 2 ) (57) 

The DFT transformation (9) is thus given by 


15 


20 


25 


30 


1 * 1 

' l 

°) 

Wd U 2 i-° = 

,0 

J 

*22- 1 


(54) 


35 


F” = Kf"'»-j= (z'+z 2 ,?!-* 2 ). 
V 2 


(58) 


where the complex conjugation is kept for consistency, and 
the numeral “1” beneath the kernel symbol in Schouten’s 
notation labels this case as a distinct solution, different from 40 
other solutions, say 


As a consistency check on both the solution for the a B of 
(56). and the legitimacy of constructing the inverse transfor- 
mation by use of the adjoint operation in (20), i.e., 


a 

2 


ir 




45 


The labeling of kernel symbols in this way is used extensively 
in the eigen-problem discussed later. 


we can check the inverse transformation calculated by way of 
the adjoint: 


ro t = = 


l l l i n l l l ri l 

i — i J V2” L i -iJvTU -i 

a Tlk W*j ° 7 ll 



( 59 ) 


55 


We comment that since (54) is positive definite, as dis- 
cussed at greater length below in the Section titled General- 
ized DFT, there will be no distinction between the other 
various quantities that may be derived from the a fi in the U 2 , 60 
and thus, no numerical distinction exists for components 
labeled by upper/lower indices, or even the complex conju- 
gation of indices. This is seen in the solution (54) above. 

As illustrated above there are an infinite number of solu- 65 
tions for the a B . Several specific solutions are shown in (55) 
based on various choices for the x 21 and x 22 components: 


in agreement with (48). The inverse DFT transformation is 
also checked for consistency giving 



--(z l ,z 2 ) = f n . 


(60) 


Several other quantity are derived from the f ' using the a R : 

/=( zV^^z'+zV-z 2 ) 


( 61 ) 
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and comparing the RHS’s in (61), as expected we have i- :: 
(f„), showing consistency using the hybrid-index notation. 
Correspondingly, various other F m quantities are given by 


F m = -}=(z' + z 2 , z 1 -z 2 ); F m = a ml F k = V2 (z 1 , z 2 ) 
F m = -j= (z T + z 2 , z T - z 2 ); F m = amn r = V2 (z T , Z 2 ) 


(62) 


and components for the W" m and the W-„ are also calculated: 


W = cf^WT 


VT 


(63) 


VT 


Wmn — 


V2 


V2 , 


Using the fj from (61) and then comparing with the F"' in (62) 
we can check the forward DFT for consistency using the W" 1 
e of (63): 


w mn f„ = 


l 

7? 

0 



V2 



(64) 


The adjoint operation on the W mn and W s „ of (63) is similarly 
verified: 


details on this below) there can in principle at least, be dif- 
ferent solutions for the eigen-problem corresponding to dif- 
ferent solutions for the a B . Note that all earlier work on the 
DFT eigen-problem uses the a B =8 B solution without consid- 
5 ering the possibility that the a B might be defined in a more 
general way based on unitary invariance. Thus, different solu- 
tions for the a B may lead to possibly different DFT eigen- 
quantities, but all the while preserving unitary invariance and 
the hennitian character of the a B . Thus, from the perspective 
10 of the transformation theory of complex vector spaces, we 
can regard some properties of the DFT as deriving from the 
characteristics of the a B , i.e., that (a) a B is invariant under 
unitary transformations and (b) is hermitian. Condition (b) 
15 really derives from (a) but it is convenient when solving for 
the a B components (as illustrated above) to regard (b) as an 
independent condition. 

Similarly as in the 2x2 example, we can solve for the a B 
corresponding to higher dimensional data sets. Of special 
20 interest are solutions in which the a B are complex valued since 
these solutions provide further consistency checks on the 
proper use of Schouten’s’ notation. We proceed by deriving 
the 4x4 W.j"' and its inverse, 

25 

-l 

rr 


from (46): 


30 


it 1 l ■ 

1 1 — / — 1 i 
21—1 1 -1 ’ 
J i -1 -i , 


( 68 ) 



such that 


W^tv, =<5” 

1 k 


Finally, Parseval’s theorem is demonstrated: 

a-J"F"=a 1 J i f= Iz 1 1 2 - Iz 2 I s , 


45 


( 66 ) 


50 


wr = w = w : 


-continued 
•11 1 1 
11 / -I —i 

2 1-1 1 -1 ’ 

, 1 -Z -1 i , 


(67) and a general 4x4 metric is defined by 


which is not in Euclidean form but is nevertheless invariant. 

Conunents on a General DFT 

Further restrictions may be placed on the solution values 
for the x 21 and x 22 in (53), in accordance with the classifica- 
tion of the a B as being positive definite, negative definite, or 
indefinite, corresponding to the number of values (or signa- 
ture)^, <0, or mixed, respectively, along the diagonal of the 
a H ). Note that when a B is definite that there is no distinction 
between the various quantities that may be derived from a B in 
the U 2 and thus, no numerical distinction exists for compo- 
nents labeled by upper/lower indices, or even the complex 
conjugation of indices; these notations are merely formal in 
this instance. We comment further that since the aj ? are 
required when solving the eigen-problem in general (more 


55 

Z T1 %2 Z T3 Z T4 
Z 21 Zj 2 2 23 Z 24 

a kl 

%1 %2 %3 %4 

. %1 %2 4ji Z44 

60 


(69) 


The condition for unitary invariance (35) then gives a simul- 
taneous system of 16 equations and 32 unknowns which can 
be further reduced to 6 free parameters by imposing the 
65 hermitian condition (37) on the a B . After considerable algebra 
the general solution to (35) for the a B , constrained by the 
hermitian condition (37) is given by 



15 
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*33 + 2(* 41 + * 43 ) *41 -iy 4 3 (*33 +*41 -*42 +*43 -*44 + 2iy 43 ) *41 - '.V43 ) 

*41 + V43 *44 *43 + ^43 *42 

(* 33 + *41 - *42 + *43 - *44 + 2iy 43 ) * 43 - fy 43 *33 *43 - V43 

*43 + ^43 *42 *43 + ^43 *44 


(70) 


As before, there are an infinite number of both real and 
complex valued 4x4 solutions for the a B in terms of the 6 free 10 
parameters in (70). Several example solutions corresponding 
to basic numerical values for the x and y components are 
listed in (71): 


15 


and then the DFT transformation (9) is given by 


F" = K,f" = ; 


z 1 + z 2 + Z 3 + z 4 
z' - iz 2 - z 3 + iz 4 
z 1 - z 2 + z 3 - z 4 


(75) 



n 

0 

0 

O’ 



0 

1 

0 

0 


X33-»U 44 -»1 - 
x 4l ->0^ 42 ->0 

0 

0 

1 

0 

’ 

X43->°’y43->° 

,0 

0 

0 

1 , 




1 

-i 

2 i 

-/' 



/ 

1 

i 

0 

2 : 33 ^ 1 ^ 44^1 = 
* 4 1 “* 0 ,x 4 2-»0 

- 

2 i 

-i 

1 

-i 

Ar 4 3 - * 0 ,y 43 -»l 


i 

0 

i 

1 , 



1 

-1 

0 

-1 ■ 


-1 

1 

1 

0 

*33->l, *44^1 = 

0 

1 

1 

1 

<41, ^-1*42^0 

i 

x 43 ~* l ’ y 43 ~*^ 

. -1 

0 

1 

1 , 


(71) 


{z'+ie-e-iz 4 } 


Applying the inverse DFT transformation shows that as 
20 expected: 


25 


W m F m = - 

4 


fl 1 1 1 ' 

' l' 4- z 2 + z 2 + Z A ' 

1 / — 1 -i 

Z l - /z 2 - Z 3 4 - fe 4 

1 -11-1 

z‘-z 2 +z 3 -z 4 

< 1 — i —1 / 

1 , • 2 3 -4 


z + iz - z - IZ 


(76) 


30 


= (z 1 , Z 2 , Z 3 , Z 4 ) 
= /"■ 



-l 

-1 

-1 

-1 v 


-1 

1 

0 

0 

£ 

"i i 
£ 1 
\L 
11 

-1 

0 

1 

0 

| x 43 -»0,y 4 3^0 | 

k -1 

0 

0 

1 , 


Other quantities are derived from the f ' using the a B of (72), 
35 e.g., consider 


40 


h ~ a tnf" ~ 2 


Zl - 1Z2 + 2 iZ3 - iZ4 ' 
!Zl + Z2 + 1Z3 
-2;'Z1 - iz 2 + Z3 - 1Z4 
!Zl + IZ3 + Z4 1 


(77) 


Once again we pick a solution and then illustrate a few of 
the basic quantities that are derived for the DFT in Schouten’ s 

notation. For generality, consider the complex-valued solu- and we further check that f ™’ ^ calculating: 
tion of (7 1 ) for the a B : 45 


a 


^ 2 Tit 


1 — / 2 i —i ' 

i 1/0 
— 2 / —i 1 —i 
i Oil. 


and then the inverse is given by 


(72) 


50 


55 


%„</*) = a fnf 
= 7n 

' Zl + iz 2 - 2i'Z3 + iZ4 
_1 -iZl+Z2-'Z3 

2 -2 ill + 1Z2 + Z3 + iZ4 

-Hi -ii} +Za 


(78) 


( 1 

-2-i 

-2 + 2/ -2- 

-2 + 2/ 

5 

2 + i 

-2 

-2 - 2 i 

2-/ 

1 

2-/ 

^ — 2 + / 

-2 

2 + i 

5 


To illustrate some specific quantities we define a general 
complex “data” vector: 


Comparing f- in (77) to the complex conjugation of F„ in (78) 
we see that as expected: 

,/HU, (79) 

demonstrating consistency using the hybrid-index notation. 
Correspondingly, we can lower indices on the F m to get 

Cf = < 8 °) 


/’=(z 1 ,z 2 ,z 3 ,z 4 ), 


(74) 
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-continued 

V +(1-20z 2 + (1+4()z 3 + (1-20z 4 

1 ( 1 + 2iV - iz 2 - ( 1 - 2i)iz 3 + i z 4 

2 (1 - 40z* - (1 + 2i)z 2 + z 3 - (1 + 2pz 4 

(1 + 2i)z‘ + tz 2 + (2i - l)z 3 - iz 4 


of a genera] rank-2 quantity, T-„. The scalar values, 


5 


are the eigenvalues and the 


Parseval’s theorem is invariant as expected: 


to 


a 

2 




F F n = 


?n/ f /' = lz 1 | 2 + |z 2 | 2 


+ lzT + lz 4 l 


■4| 2 


( 81 ) 


but the Euclidean form in (81 ) deriving from this particular 


are the eigenvectors. The “a” index below the kernel object, in 
this case 


V 

a 


a 

2 


mn 


20 and 


is not necessarily expected. As pointed out earlier, Parseval’s 
theorem can take different forms for the different a B solutions. a- 

E.g., considering solution 4 of the a s/ in (7 1 ) we have 0 


V" 




does not refer to a component, but rather, the a-th eigenvector, 

( 82 ) 


= -\z l \ 2 + \z 2 \ 2 + k 3 | 2 + \z 4 \ 2 - 2 z l (z 2 + z 3 + z 4 ), 


30 


which is invariant but not real valued. 

Eigenvalue-Eigenvector Problem 

Given the formalism above, we can examine the eigen- 35 
value-eigenvector (eigen-problem) from a more general 
viewpoint than what is normally discussed. The additional 
generality stems from the possibility that the metric is not 
necessarily the Kronecker delta, constrained only by its 4(| 
invariance under unitary transformations and hennitian sym- 
metry. We have shown from a general viewpoint that the a 
e= 8 n solution is an important “base” solution without addi- 
tional worries that perhaps something was overlooked. We 
have also shown that other solutions for the a I; exist based 4 _ 
simply on unitary invariance, and thus, more general calcu- 
lations of the DFT and its various quantities are easily 
extended to arbitrary basis sets. Some advantages in calcula- 
tion may be possible in these other general basis sets and we 
will reserve this topic for a future set of notes. For now, we 5Q 
examine a specific basis set that follows by solution to the 
eigen-problem as discussed below. 

Before examining the specifics of the eigen-problem for 
the DFT, some generalities are discussed. The eigenvalue- 
eigenvector problem seeks solutions to the following set 55 
equations: 

Tmnf = a-Vm, ( 83 ) 

60 

for the vectors 


corresponding to the a-th eigenvalue, 


IT 

a 

as separate and distinct solutions. We therefore choose indi- 
ces from the set, {a, b, c, d}, to label the kernel objects 
corresponding to distinct eigen-solutions, and then reserve 
the index range, {k, 1 , m, n}, to label components of a given 
basis. 

Continuing to the DFT application we substitute W Sm of 
(26) into (83): 


W mn v n =o-v m , (84) 


and make special note that we have switched to a common 
kernel symbol for the Fourier transform pair, (F m , f)-*(v;j, 
v"), with the substitutions: 

^/■andfipv^. (85) 

Looking ahead, we mention now that the solution to (84) 
results in a transformation of the W-„ into a diagonal basis by 
use of the similarity transform. In the context of the present 
application, and noting expression (9) for the Fourier trans- 
form, (84) then has an interesting geometrical interpretation 
as solving for the Fourier transform pair, (v-, v”), that are 
“parallel” to one another, or equivalently, left invariant under 
the action of the DFT (aside from the scaling factor cr). To give 
greater emphasis to this last comment, consider the overall 
structure of the DFT transformation in a diagonal basis for the 
4x4 case (for now using the notation of (9) to emphasize this 
last point): 
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W mn v n = crv m 


o- 0 0 0 

l 

0 o- 0 0 

2 

0 0 o- 0 

3 

0 0 0 O- 


r/ if ' 

/ : 

/ : 

/■ 


( (cr /cr)/ 1 ' 

| \ 1 / a / 

I a ) 

/ a/ 

I (cr / <r\f 4 ' 
V \4 / flj 


20 


review some of these basic results before examining the 
(86) diagonal DFT calculation in detail. Since 


w*=l 


(91) 


the eigenvalues for the DFT matrix for N>3 are non-unique 
and are multiples of the so-called “roots of unity:” 


15 


V 

? / 2 ' 

P 2 ' 

T/ 2 ' 

Fi> 

Ft 
Fy 
F« 

where the primes on indices denote the data and its transform 20 
in the diagonal basis as discussed further below. Thus as 
mentioned above, the Fourier transform pair are co-linear 
when calculated in the diagonal basis: 


O- = e"“ /2 ; a = 1 4; => cr e {±1, ±(). 


Specifically, (92) then becomes 


(92) 


cr — 1, cr — — 1, cr = / cr = - 

12 3 4 


(93) 


The eigenvalues of the DFT transform are not all distinct 
for N>4, so there are many possible sets of eigenvectors 
satisfying (89). As a consequence of having these non-unique 
eigenvalues, the corresponding eigenvectors for these degen- 
erate eigenvalues are not necessarily orthogonal. The impor- 
tance of having an orthogonal set of eigenvectors for the 
(87) 25 present application is not known, but should this turn out to be 
an issue, the Gramm-Schmidt procedure can always be 
applied or an orthogonal set can be obtained. Further conse- 
quences and examples of this analysis for higher dimensional 
DFTs will be discussed in a future set of notes. For now we 
examine a low-dimensional example to map out some of the 
basic features. 


30 


differing only by a scalar factor. This example illustrates the 
geometrical significance of the DFT transform in a diagonal 
basis and motivates using the same kernel symbol for the 
Fourier transform pair in (85), since the quantities are co- 
linear. This difference in notation is to be contrasted with the 
approach taken so far where the kernel symbols have been 
kept distinct, and the W.„ m transformation was regarded as a 
point transformation. 

Continuing with the eigen-problem calculation we use the 
metric tensor to show that (84) is equivalently expressed in 
the form: 


35 


( Wmn CQmnjV — 0 , 


( 88 ) 


Eigenvalues 

Considering M=N=4, the DFT transform and its inverse is 
given in (68) which is reproduced below: 


40 


W% = - 

k 2 


1111 
1 -i -1 i 
1-11-1 
1 i -1 -i 


45 


: W % m k 


1111 
1 i -1 -i 
1-11-1 
1 -i -1 i 


or also as 


Considering solutions of (90) we have 


50 


W™ - crd™ v" = 0. 


(89) 


del( w^ - tro™ j = 0, 


And thus the eigen-problem, when cast in the form of (89), 
does not depend on a specific choice for the metric. Thus we 
solve the DFT eigen-problem in the form of (89). The eigen- 
values for the W.„ m are then solutions to 


and for the 4x4 example has the solution: 


(T - 1 )% + 1 )(o + ‘) = ° => {T= 1 ’T = -1, T =1, T = - ‘}’ 


( 94 ) 


det W"‘ - (rd™ 


(90) 


with one degenerate eigenvalue, 


Eigen-Problem for the DFT Matrix 65 

The eigenvector-eigenvalue problem for the DFT matrix in 
the form of (89) has been examined extensively and so we 


- = 1 . 
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Multiplicity of Eigenvalues for the DFT Matrix 
As discussed above, the significance of (91) for the DFT 
eigenvalue problem is that degenerate eigenvalues appear for 
N>4. The “multiplicities” of these eigenvalues are defined as 
the number of times each eigenvalue appears. For the 4x4 
example in (94) we therefore have 

m 1 =2, m 2 -=l, m 3 =l, m 4 =0, (95) 

where the m k ordering corresponds to the ordering of (93). 
These results can be generalized to higher dimensions by the 
following Table of multiplicities for each eigenvalue: 

TABLE 1 


which simplifies to 


5 


to 


-v 1 + v 2 + v 3 + v 4 
llil 


v 1 -(2 + /)v 2 -v 3 iv 4 

l ill 



v l + iv 2 - v 3 - (2 + 1 > 4 
ill l 


The solution to (100) for 


(°'l 

0 

0 

. 0 , 


(100) 


M = N 

m i 

for cr = 1: 

l 

m 2 

for cr = - 1 : 
2 

m 3 

for cr = i: 

3 

m 4 

for cr = — i: 

4 

4m 

m + 1 

m 

m 

m- 1 

4m + 1 

m + 1 

m 

m 

m 

4m + 2 

m + 1 

m + 1 

m 

m 

4m + 3 

m + 1 

m + 1 

m + 1 

m 


By inspection of Table 1 , the 4x4 case corresponds to the m= 1 
value in the top row giving the multiplicities of eigenvalues 
listed in (95). Similarly, it can be shown that for the 5x5 case, 
the eigenvalues are solutions of 


15 


20 


25 


v" 

l 


is therefore given in terms of two arbitrary components, say 


v 3 and v 4 : 

i i 


(v 1 , V 2 , V 3 , V 4 ) = (v 3 + 2v 4 , I' 4 , V 3 , V 4 V 
' i i l i / Vi 1111/ 


( 101 ) 


(cr-lffo-H 


1 <7**1 


( 96 ) 


> (cr = 1, cr = - 
1 1 2 
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We point out that two components are required because of 
the two-fold degeneracy of the 


corresponding to the 2 nd row of the Table with m=l. From 
(96) (or equivalently Table 1) the multiplicities for the 5x5 35 
case are thus 

w 1 =2, m 2 =l, m 3 =l, m 4 =l. (97) 

To illustrate a higher dimensional example, consider say 
M=N=131=4.32+3, corresponding to the third row of the 
Table giving the following multiplicities: 

w t = 33, m 2 = 22, m 3 = 22, m 4 =32. (98) 

Eigenvectors 

Given the eigenvalues in (94) the eigenvectors correspond- 
ing to each eigenvalue are constructed using (89) . It is instrac- 4 - 
five to consider this process for the eigenvectors correspond- 
ing to the degenerate first and third eigenvalues, 


1 3 

The first eigenvector is obtained as a solution to the set of 
linear equations: 

(99) 


( 1 - cr 

l 

1 1 1 ' 

(v 1 

1 


(0\ 

1 

—i — cr -1 i 

l 

V 2 

1 


0 

1 

-1 l-o- -1 

V 3 

- 

0 

1 

i —1 —i — cr 

l 

1 

v 4 

1 


.0, 


cr = cr = 1 
1 3 

eigenvalue. In general, the eigenvectors corresponding to 
unique eigenvalues are specified in terms of a single compo- 
nent. Continuing we arbitrarily choose 



and 


v 3 =0, 

1 

to obtain a realization for the first eigenvector: 



= (2, 1, 0, 1). 


Note also that since, 


cr = cr = 1, 

1 3 

are degenerate, the solution to (99) is trivially the same as 
(100). Therefore, a solution for 
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of (71): 


5 

follows immediately by choosing complementary but arbi- 
trary numerical values for the 


'3 0 1 

/II j. . tt . 020 

(d„ v ) = ^f„ v=,5 »f„ v = 1 0 1 

v 0 0 0 


0\ 

0 

0 ’ 


u 


(105) 


10 where the off-diagonal entries correspond to the 


and 


15 


cr = cr = 1 

1 3 


value. For reference, we can still construct an orthonormal set 
from (104): 


in (102), e.g., 


v 4 = 0 


3 


20 


25 


■2' 


-t ■ 


1 ' 


f 0 ’ 

1 

2 2 

1 

3 2\fJ 

-1 

. 1 

-l 

0 

1 

3 

, U — — 

4 VT 

0 

. 1 , 


1 , 


, -1 , 


, l , 


(106) 


v 3 = 1: 

3 


v 1 , V 2 , v 3 , V 4 ) = (v 3 + 2v 4 , v 4 , v 3 , v 4 } v" = (1, 0, 1, 0) 

.3 3 3 3 / V3 3 3 3 3 / 3 


and verily as expected that 


( 103 ) 


30 = 6 mn . 


007) 


Similarly, we can solve for the remaining eigenvectors 
corresponding to their eigenvalues, 


cr = - 1 

2 


Similarity Transform 

35 Given the eigenvectors of (104) we form the “similarity” 
transform using the eigenvectors as columns. We denote this 
as a transformation of basis and we find by construction that 
the similarity transform, S m m ', and its inverse, 


and 


40 


45 


are given by 


In summary, a set of eigenvectors (but not the only set) for the 
4x4 DFT transformation is summarized below: 
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1 
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0 
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'1' 


' 0 ' 

1 

, v” = 

1 

, v n = 

0 

> v” = 

-1 

0 

2 

1 

3 

1 

4 

0 

4, 


< 1 , 


c0. 


5 1 , 


Since two of the eigenvalues are degenerate, (102) and 
(104) collectively do not form an orthogonal system as can be 60 
verified by forming the inner product between any two vec- 
tors using the “base” solution, 


When using the kernel -index method for transformations 
of basis, say x m =A m m x m , that the A m m ' quantities do not use 
dots above or below an index to mark index placement. This 
differs from the notation for point transformations, which 
does use the dots e.g., F m =W.„ m f'. The reason for this is to 
discourage the “raising” or “lowering” of indices on the A m m 
to avoid an inconsistency when using the metric tensor. This 
inconsistency can occur because the metric a B also transforms 
under the A m m ', and thus has different components in the new 
coordinate system. For point transformations however, no 
coordinate transformation occurs and thus indices can be 
raised or lowered using a single (unambiguous) metric. 

We then verify that applying the similarity transformation 
to the W.„ m gives as expected: 



25 


US 9,015,149 B2 


26 


sws~ l 


= sz w™s „ 


(O- 0 0 0 

1 

0 O- 0 0 

2 

0 0 o- 0 

3 

0 0 0 £7- 

4 


1 0 0 0 ' 
0-100 
0 0 10 
0 0 0 -i. 


(109) _! , 

S n . 

5 + * 
to give 


As a check on index ordering we run the following line of 
Mathematica code with Si defined as the inverse of the Sin 
(108): 

TsLble[Sum[S[[m,ump]JW[[m,nJJSi[/lnp,n]],{n,l,M}, 

{ump,l,M}, {Inp, l,M}] (110) 

where ump denotes “upper-m-prime” and lnp is “lower-m- 
prime.” 

It is of further interest to point out that if the similarity 
transform is constructed using the orthonormal eigenvectors 
of (106) as columns (and correspondingly it inverse), then the 0(| 
application of the “orthononnal” similarity transform to the 
W.,," 1 also gives (109). Therefore, the similarity transform that 
leaves the W.„ m in diagonal form is certainly not unique, so 
many possible choices can be made for the S,," 1 ’. So as indi- 
cated by calculations thus far, it does not seem to matter 
whether or not the eigenvectors forming the S„ m ’ are orthogo- 
nal, and this can greatly simplify the initial calculations 
needed for calculating in the diagonal basis. 

DFT Calculation in Diagonal Form 

The diagonal form of the DFT in (109) shows that the DFT , f) 
calculation can be performed both efficiently and elegantly as 
simply a multiplication of the entries on the main diagonal 
(the eigenvalues), to the corresponding data components that 
have been pre-multiplied using the similarity matrix. This 
means that after pre-multiplying the data, the DFT transform , . 
is reduced to a vector product consisting of N operations 
rather than a matrix multiply with N 2 operations (or N log N 
when using the FFT). In reality though, this reduction in 
computation cost is only partly true in terms of the perfor- 
mance gains, i.e., the diagonal DFT representation does obvi- 
ously lead to a speedy calculation for one part of the DFT 
operation, but in general, this gain is offset by the added 
computational cost of having to pre-transform the data and 
then back-transform to bring the Fourier transform into the 
original basis as illustrated below. But the key is that these 45 
latter operations may not need to be performed at every step 
in an iterative transform application, such as with phase 
retrieval. Additional details on this calculation approach are 
outlined below. 

To begin, we write down the DFT calculation in the diago- 5Q 
nal basis after switching the notation back to distinct kernel 
symbols for the data and its transform: 

r'=W. n m 'f, (ill) 

where the m' represents the transformation of basis under the 55 
action of the similarity transform as illustrated above in ( 1 09). 
Since the only non-zero values for the W are along the 
diagonal, we form a vector, 


/" = S’ /». 


(112) 


The diagonal DFT is then the product of the N components 
of the 


and the f' defined by (112): 



(113) 


where no sum is implied over the (m',m') pair. To complete the 
calculation, the “diagonal” DFT components in (1 13) need to 
be transformed back to the original DFT basis, F m . To obtain 
this result we express the LF1S of (1 1 1) in terms of compo- 
nents in the original basis: 


F m = s™ F m , 


(114) 


and therefore the original Fourier transform components are 
readily obtained by back -transforming the F” 1 ' of (114) using 
the similarity transform: 


F" = 5", /." _ 5» pm - fit pm 


(115) 


The inverse DFT is also calculated in a diagonal basis since 
the inverse of a diagonal matrix is found by populating its 
diagonal entries with the inverse of the individual compo- 
nents. Therefore, given the the inverse is simply (no 

sum is implied on the RHS over the (in', m') pair): 



(116) 


and then inverse DFT transformation in the diagonal basis 
gives as expected: 


60 = ('h\z r ' 


(117) 


form the eigenvalues (i.e., the diagonal entries of the W.„, m ) 
as in (94). Continuing with the calculation of (1 1 1, the start- 
ing data, f fo are pre-multiplied using the inverse similarity 
transform, 


The overall point is that for applications where it is feasible to 
perform the DFT in a diagonal basis, the similarity transform 
and its inverse can be pre-calculated. The multiplication of 
the similarity transform with the data can also be pre-calcu- 
lated. Then at an appropriate point, the DFT results can be 
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transformed back to the original basis using the inverse simi- 
larity transform, preferably at the end of an iterative loop. The 
approach is su mm arized in FIG. 2. 

Example 5 

4x4 Diagonal DFT Calculation 

The steps above are illustrated for the 4x4 example. Con- 
sider the general complex data vector of (74): io 


/=(z 1 ,z 2 ,z 3 ,z 4 ). 


(118) 


The inverse similarity transform of (108) times the data 
vector (1 18) is given by 


15 


-i / 

--sir 


( Zl + Z2 -Z3 + Z4 
— Zl + Z2 + Z3 + Z4 
Zl - Z2 + 3z 3 - Z4 
2(-Z2 +Z4) 


(119) 


20 


and then from (94): 

25 

<r = (1,-1, 1,-0. (120) 


The DFT in the transformed basis follows by multiplica- 
tion of the eigenvalues in (120) with the f"' components in 30 
(1 19) to give the Fourier transform in the diagonal basis: 


( Zl + Z2 - Z3 + Z4 ) 


F m> =cr-f m = - 

m' 4 


Zl Z2 Z$ Za 
Zl -Z2 +3 z 3 -Z4 
-2i'(-Z 2 +Z4) 


( 121 ) 

35 


Finally, we back-transform the F m ' components in ( 1 2 1 ) by 40 
multiplying with the similarity transformation of (108) to 
give the DFT result in the original non-diagonal basis: 


F m = 5 ”, F m ' 


(122) 45 
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0 ' 


’ Zl +Z2 -Z3 +Z4 ’ 

1 

1 

0 

-1 

1 

Zl - Z2 - Z3 - Z4 

0 

1 

1 

0 

4 

Zl -Z2 +3a -Z4 

, 1 

1 

0 

1 , 


. -24-Z2+Z4) . 


tational cost of setting up the data in a transformed basis as 
illustrated above. But even with this in mind, certain applica- 
tions spend most of the time in the forward and inverse trans- 
form part of the calculations (see FIG. 2), so the main point is 
to explore applications where the initial and final data trans- 
formations are needed only once, and then subsequent DFT 
calculations can proceed as simply N multiplications of the 
entries along the main diagonal. The performance gain for a 
single DFT calculation can then be approximately log(N) 
times faster. Assuming for simplicity a power of 2 for the data 
size (N=2"), we can estimate a performance improvement for 
large n to be n log(2)~10 1 , or about an order of magnitude 
faster than the equivalent FFT calculation, and for arbitrary 
frequency spacing. 

FIG. 3 shows an example embodiment of the present inven- 
tion and is thus directed to an image-based wavefront sensing 
method of and means for focusing by determining phase 
difference between received images, removing that phase 
difference by calculating a DFT for each signal and iteratively 
multiplying the resultant until a convergence is achieved by 
accepting point-source stellar objects to recover embedded 
optical phase information, receiving data from a sensor array 
wherein the data includes embedded phase information, cal- 
culating a similarity transform and pre-transforming the data 
once at the beginning of the calculation. 

FIG. 3 further shows, in an example embodiment of the 
present invention, an image-based wavefront sensing method 
of and means for iteratively processing the data using arbi- 
trary frequency spacing in a diagonal fashion until a wave- 
front convergence is achieved to determine a phase difference 
by implementing the DFT in a 1 dimensional linear array 
using just N operations, where 0<N<Nxlog(N)<oo, wherein 
Nxlog(N) is the number of operations of an FFT counterpart, 
back-transfomiing the data at the end of the calculation in a 
single step and focusing the received images based on the 
determined phase. 

While the invention is described with reference to an exem- 
plary embodiment, it will be understood by those skilled in 
the art that various changes may be made and equivalence 
may be substituted for elements thereof without departing 
from the scope of the invention. In addition, many modifica- 
tions may be made to the teachings of the invention to adapt 
to a particular situation without departing from the scope 
thereof. Therefore, it is intended that the invention not be 
limited the embodiments disclosed for carrying out this 
invention, but that the invention includes all embodiments 
falling with the scope of the appended claims. Moreover, the 
use of the terms first, second, etc . does not denote any order of 
importance, but rather the terms first, second, etc. are used to 
distinguish one element from another. 


in agreement with the original DFT transform, F m W.„ m f ' of 
(75). Similarly, the inverse transformation can be calculated 
in the diagonal basis as discussed above using (117). 60 

The diagonal form of the DFT in (113) shows that the 
calculation can be performed, efficiently and elegantly as 
simply a vector product. This means that the DFT can be 
reduced to N multiplications of the data with the eigenvalues. 
But as emphasized above, this observation is only partly true, 65 
i.e., although the performance gains from the diagonal repre- 
sentation are indeed realized, these are offset by the compu- 


What is claimed is: 

1 . An image-based wavefront sensing method in a space 
based interstellar observatory with mirror segments of focus- 
ing by determining phase difference between received 
images, removing that phase difference by calculating a dis- 
crete Fourier transform (DFT) to produce a series of result- 
ants for a plurality of signals corresponding to received 
images and iteratively multiplying the resultants with previ- 
ous iterative resultants of each signal until a convergence is 
achieved comprising the steps of: 

providing a space based interstellar observatory with mir- 
ror segments; 

accepting point-source images to recover embedded phase 
information by said space based interstellar observatory; 
receiving data from an external sensor array wherein the 
data includes embedded phase information; 
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pre-calculating a similarity transform based on orthonor- 
mal eigenvectors as columns; 

pre-calculating a multiplication of the pre-calculated simi- 
larity transform with the data to minimize calculation 
time and processor hardware; 

iteratively processing the data in a diagonal fashion to 
determine a phase difference by implementing the DFT 
in a 1 dimensional linear array using just N operations, 
where 0<N<Nxlog(N)«», wherein Nxlog(N) is the 
number of operations of an FFT counterpart; 

back-transforming the data at the end of the calculation in 
a single step; 

minimizing pretransforming and back-transforming data 
operations; 

focusing the received images based on the determined 
phase difference; and 

iteratively realigning the mirror segments to minimize 
focusing error effects. 


30 

2. The wavefront sensing method of claim 1 wherein cal- 
culating a similarity transform further includes pre-trans- 
forming the data once at the beginning of the calculation. 

3. The wavefront sensing method of claim 2 wherein itera- 
5 tively processing further includes a single step of back-trans- 

forming the data at the end of the calculation. 

4. The wavefront sensing method of claim 3 wherein said 
point source images include stellar objects. 

5. The wavefront sensing method of claim 4 wherein said 
embedded phase information is of an optical nature. 

6. The wavefront sensing method of claim 5 wherein said 
external source includes a sensor array. 

7. The wavefront sensing method of claim 6 wherein the 
step of iteratively processing further includes the step of 
determining the phase difference by utilizing a discrete Fou- 

15 rier transform (DFT) . 

8. The wavefront sensing method of claim 7 including 
iteratively processing using arbitrary frequency spacing. 



